clear; clc; close all;

% 参数
a = 0.1;
b = 0.1;
c = 0.2;     % Fig.10(b) 中 c=0.2
k = 1.85;

% 网格范围与分辨率，可根据需要调整
xrange = linspace(-4, 4, 401);  
qrange = linspace(-4, 4, 401);

% 迭代与判定参数
maxIter = 600;     % 总迭代步数
transient = 200;   % 暂态步数
threshold = 1e3;   % 发散阈值
var_epsilon = 1e-8;% 判定定点吸引子的方差阈值

% 建立分类矩阵 colorMatrix
% 0=发散(白), 1=定点(黑), 2=混沌/周期(青)
colorMatrix = zeros(length(qrange), length(xrange));

for ix = 1:length(xrange)
    for iq = 1:length(qrange)
        x0 = xrange(ix);
        q0 = qrange(iq);
        isDiverge = false;

        % 先迭代 transient 步去除暂态
        for t = 1:transient
            [x0, q0] = mptm(x0, q0, a, b, c, k);
            if abs(x0) > threshold || abs(q0) > threshold
                isDiverge = true;
                break;
            end
        end
        
        if isDiverge
            colorMatrix(iq, ix) = 0; % 发散
            continue;
        end
        
        % 再记录剩余 maxIter-transient 步
        remainSteps = maxIter - transient;
        orbit_x = zeros(1, remainSteps);
        orbit_q = zeros(1, remainSteps);
        orbit_x(1) = x0;
        orbit_q(1) = q0;
        for t = 2:remainSteps
            [orbit_x(t), orbit_q(t)] = mptm(orbit_x(t-1), orbit_q(t-1), a, b, c, k);
            if abs(orbit_x(t)) > threshold || abs(orbit_q(t)) > threshold
                isDiverge = true;
                break;
            end
        end
        
        if isDiverge
            colorMatrix(iq, ix) = 0; % 发散
        else
            % 判断定点还是混沌/周期
            vx = var(orbit_x);
            vq = var(orbit_q);
            if vx < var_epsilon && vq < var_epsilon
                colorMatrix(iq, ix) = 1; % 定点
            else
                colorMatrix(iq, ix) = 2; % 混沌/周期
            end
        end
    end
end

figure('Name','Fig.10(b)');
hold on;
xlabel('x_0'); ylabel('q_0');
title('The basin of attraction (c=0.2, k=1.85)');

% 将分类矩阵映射到 RGB
%   0=发散(白), 1=定点(黑), 2=混沌(青)
RGB = ones(length(qrange), length(xrange), 3);  % 默认全白
for iq = 1:length(qrange)
    for ix = 1:length(xrange)
        if colorMatrix(iq, ix) == 1
            RGB(iq, ix, :) = [0, 0, 0];  % 黑
        elseif colorMatrix(iq, ix) == 2
            RGB(iq, ix, :) = [0, 1, 1];  % 青
        end
    end
end

imagesc(xrange, qrange, RGB);
